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Abstract 

Many mathematical models involve input parameters, which are 
not precisely known. Global sensitivity analysis aims to identify the 
parameters whose uncertainty has the largest impact on the variability 
of a quantity of interest (output of the model) . One of the statistical 
tools used to quantify the influence of each input variable on the output 
is the Sobol sensitivity index, which can be estimated using a large 
sample of evaluations of the output. We propose a variance reduction 
technique, based on the availability of a fast approximation of the 
output, which can enable significant computational savings when the 
output is costly to evaluate. 
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Introduction 

Many mathematical models encountered in applied sciences involve a large 
number of poorly-known parameters as inputs. It is important for the prac- 
titioner to assess the impact of this uncertainty on the model output. An 
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aspect of this assessment is sensitivity analysis, which aims to identify the 
most sensitive parameters, that is, parameters having the largest influence 
of the output. In global stochastic sensitivity analysis (see for example jl| 
and references therein) the input variables are assumed to be independent 
random variables. Their probability distributions account for the practi- 
tioner's belief about the input uncertainty. This turns the model output 
into a random variable, whose total variance can be split down into different 
partial variances (this is the so-called Hoeffding decomposition, see 
Each of these partial variances measures the uncertainty on the output in- 
duced by each input variable uncertainty. By considering the ratio of each 
partial variance to the total variance, we obtain a measure of importance 
for each input variable that is called the Sohol index or sensitivity index of 
the variable ll|]; the most sensitive parameters can then be identified and 
ranked as the parameters with the largest Sobol indices. 
Once the Sobol indices have been defined, the question of their effective 
computation or estimation remains open. In practice, one has to estimate (in 
a statistical sense) those indices using a finite sample (of size typically in the 
order of hundreds of thousands) of evaluations of model outputs [1]. Indeed, 
many Monte Carlo or quasi Monte Carlo approaches have been developed by 
the experimental sciences and engineering communities. Such an approach 



is the Sobol pick- freeze (SPF) scheme (see [111. Il2l]). In SPF a Sobol index is 
viewed as the regression coefficient between the output of the model and its 
pick-freezed replication. This replication is obtained by holding the value of 
the variable of interest (frozen variable) and by sampling the other variables 
(picked variables). The sampled replications are then combined to produce 
an estimator of the Sobol index. 

The SPF method requires many (typically, around one thousand times the 
number of input variables) evaluations of the model output. In many in- 
teresting cases, an evaluation of the model output is made by a complex 
computer code (for instance, a numerical partial differential equation solv- 
ing algorithm) whose running time is not negligible (typically in the order 
of a second or a minute) for one single evaluation. When thousands of such 
evaluations have to be made, one generally replaces the original exact model 
by a faster-to-run metamodel (also known in the literature as surrogate model 
or response surface [1]) which is an approximation of the true model. Well- 



known metamodels include Kriging [10], polynomial chaos expansion 13|] 
and reduced bases 0) IB] > to name a few. From a multifidelity point of view, 
the metamodel can also be viewed as a "coarse" (low-fidelity) version of the 
code; the metamodel, seen as a coarse version, may also be a "degraded" 
version of the code: for instance, it may be a solver for a simplified model 
(either mathematically simplified, or discretized on a coarser grid), an inte- 
grator for a function of lesser precision, or an optimizer stopped before its 
full convergence. In this paper, we designate by "coarse approximation" any 
of the above approximations (metamodels and degraded versions). When 
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using a coarse approximation for sensitivity analysis, the original model is 
generally used only to define the metamodel, and not to perform the Sobol 
index estimation. This leads to a necessity of measuring the difference be- 
tween the model and its approximation in order to certify the sensitivity 
index estimation [1,0]. To our best knowledge, no approach for using both 
metamodel and model evaluations to estimate Sobol indices have been pro- 
posed yet. 

In this work, we propose an approach, based on the asymptotic properties 
of the SPF scheme studied in [6] to optimally combine evaluations of the 
original model and evaluations of its approximation, in order to produce an 
asymptotically-justified confidence interval for the Sobol index of the original 
model. Our approach is inspired by the quasi-control variate method j3| 
which has been developed for Monte-Carlo estimation of means. 
This paper is organized as follows: in the first section, we begin by setting 
up the notations and the context of the paper. Then we define the Sobol 
index estimator we wish to study. The main result is Theorem 11.11 which 
provides an asymptotic method to estimate a confidence interval for a Sobol 
index. The second section is a numerical illustration on a particular (but 
representative) kind of model output. 



1 Motivation and definition of the estimator 
1.1 Notation and context 

We begin by setting up the usual notations in the sensitivity analysis con- 
texts. The output of interest is a random variable Y, which is a deterministic 
function r/ : M of the random inputs X £ and Z eW^: 

Y = r]{X,Z), 

where pi and p2 are integers, and p = pi + P2- 

We assume that X and Z are independent random variables and that Y has 
a finite and nonzero variance. We are interested in the (closed) Sobol index 
with respect to X, defined by: 

^ Var(E(y|X)) 
Vary 

This index, which is between and 1, quantifies the influence of the X input 
on the output Y: a value of S that is close to 1 indicates that X is highly 
influential on Y. 

The pick-freeze method [()] expresses S using a covariance: 
where Z' is an independent copy of Z. 



3 



This expression leads to different Monte-Carlo estimators of S. For instance, 
the following estimator is studied in 
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where, (yi)i=i,...,Af and (l'/)j=i,...,Ar are independent samples of F (resp. Y'), 
and, as in the rest of the paper, all sums are for i from 1 to A^. 
It is shown [op.cit.. Proposition 2.2] that (T/v)Ar is asymptotically normal, 
with variance a^^/N, where: 

2 _ Var ((y - E{Y)){Y' - E{Y)) - S/2 {{Y - E{Y) f + {Y' - E{Y) f)) 
'^^''^ ~ (Var(y))2 ' 

(1) 

and [op.cit.. Proposition 2.5] that this asymptotic variance is minimal among 
regular estimators that are functions of realizations of exchangeable {Y, Y') 
pairs. 

Note that a realization of the estimator, for a finite sample size A^, can 
be computed by making 2N evaluations of the t] function. 
In this paper, we suppose that we can evaluate, in addition to the rj function, 
an approximation r]c : W ^ M oi the r/ function (the c index is for coarse). 
The usage of such an approximation has been motivated in the Introduction. 
A concrete and ubiquitous example of r] and rjc will be presented in the 
next section. In the following section, we motivate and study our variance- 
reduced estimator of S. 



1.2 Variance-reduced estimator 

Let: 

Y, = r],iX,Z), Y^ = 7],{X,Z'), 

and {Yci)i=i^...^N, iYci)i=i,...,N be A^-samples of (resp. Y^). 
The estimator: 



Tn 



N ^ ^^^-^ ci \ N /-^ 2 
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consistently estimates the Sobol index of the coarse model: 

Var{E{Y,\X)) 



Vary, 



by using 2N evaluations of rjc- 
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As mentioned in the introduction, our objective is to combine evaluations 

of T] and 77c to estimate S" at a smaller cost than an estimation that would 

be performed from evaluations of ij only. 

We take a function ^ : N — )• N. 

It is clear that the estimator Ej\f defined by: 

consistently estimates E = S — Sc, and that a realization of E^- can be 
obtained using 2ip{N) evaluations of rjc and 2ip{N) evaluations of ry. 
We propose a natural estimator of S based on Tjy and i?Ar, inspired by the 
quasi-control variate method is thus: 

Vn = Tn + En- 

This estimator can be computed by making 2(A^ + ^lJ{N)) evaluations of r]c 
and 2ip{N) evaluations of r/. As an evaluation of r] is more costly than one of 
r/c, one can expect a computational gain if ip{N) < N, and if the asymptotic 
variance of (Vat) is less than the asymptotic variance of (T^^^^), so that 
asymptotic confidence intervals built upon Vn are more precise than those 
built on r^(jv) ^^one. 

The following theorem gives a method for estimating (conservative) asymp- 
totic confidence intervals using Vn. We denote by $ the cumulative distri- 
bution function of the Gaussian with zero mean and unit variance, and by 
its inverse. 

Theorem 1.1. Suppose that lim^v^+cxi = +00. 
Then, for any and Uc in ]0, 1[; 

lim P |5 - Va^I < qioLe) Zl—= + g(ac)4= > 1 - ("e + "c), 

for: 

/X ^1/ , ^ 2 Var(Ac-Sc/2) 

2 _ 2 Var(A-P/2) 2Cov(^, A^) - (Cov(A, B^) + CovjB, A^)) + CovjB, Bc)/2 
(Vary)2 Vary VarY^ ' 

where A, B, Ac, Be are the following random variables: 

A = {Y - E(y))(y' - E(y)), b = s[{y - w.{Y)f + (y' - ^(y))^] , 
Ac = {Yc - E(ye))(y; - E(ye)), Bc = Sc [{Yc - E(yc))2 + (y/ - ^{Yc)f] . 

The same holds when cTc and are replaced by any consistent estimators. 

Sketch of proof . Follow the proof of fil]. Proposition 2.2, and apply the 6- 
method to (T/v,T'^^^p to get the asymptotic variance of (E^)- 
Then use that for any ei, 62, 

{\Vn -S\>ei + 62} C {|r^ - Sc\ > ei} U {\En - E\ > 62} . □ 



5 



1.3 Choice of V', cte and ac 



To convert the theorem above into a practical procedure, it remains to choose 
the parameters "0; Oe and so as to minimize the overah computational 
time. 

We will assume that one evaluation of ?] as a unit cost, and that an evaluation 
of rjc has cost < p < 1. We also set iIj{N) = [/iA^] , where // g]0, 1[ is to be 
found, and [•] is the "ceiling" function. 

We choose a target risk level a g]0, 1[ and a target length L for the confidence 
interval of Theorem 11.11 

It is clear these constraints force ac in function of and a: 

Oc = a*{ae) = 1 — (a + ae), 
and that N has to satisfy: 

N > iV*(ae,/i) := (^^^ + g(«*(«e))^c)' 

We approximate 'tp{N*) by fj,N* . The cost of the required evaluations of rj 
and r]c is thus, in the general case: 

Cost(ae) fJ-) = 2A^*(ae5 Ai) (2/^ + p) , 

corresponding to the ip{N) evaluations of r] and the il^{N) + N evaluations 
of rjc- 

However, in some settings, the computations made to compute r]c can be 
reused to compute rj, allowing to evaluate rjc and rj for a unit cost, leading 
to: 

CostHier(ae, /") = 2N*{ae, jj) {fi + p) . 

Such a "hierarchical" property is beneficial to our estimation scheme and 
occurs naturally for some ij, as we will see in the numerical illustration 
section. 

Now, one would obviously choose and p so as to minimize the cost 
Cost(ae> p) (or, depending on the case at hand, CostHier(ae) p))- In practice, 
this is not possible, as cJe and Uc are unknown. Hence, approximately opti- 
mal parameters are found by empirically estimating these quantities, based 
on a small sample of realizations of y, Y' , Yc and Y^. This gives rise to S* 
and p*, and an estimated optimal costs: 

Cost («*,/!*) and CostHier(S*, /I*). 

2 Numerical illustration 

2.1 Model set-up 

In financial mathematics, the Heston model [3] is the following stochastic 
differential model for the price of a risky asset {St)t>o as function of the 
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Name 


Interpretation 


Min. 


Max. 




Initial volatility 


.2 


.25 


K 


Volatility convergence rate 





3 


9 


Volatility limit 


.2 


.22 


r 


Correlation between Brownians 


-1 


1 




Volatility of the volatility 





.4 


R 


Risk-free rate 


.08 


1.1 



Table 1: Distributions and interpretations of the input parameters, 
time t > 0: 

( dSt = fiSt dt + ^tSt dW/ 

\ dut = K{e-iyt)dt + ^^tdWt^ ' 

where (TV/)t>o and {W^)t>o are standard Brownian motions (under the 
risk-neutral probability measure Q) whose correlation is r S [0, 1]. 
We are interested in the price of an European call option of maturity T and 
strike K, which is given by e~^'^KQ [{St — 

Although a semi-analytical formula is available for the fast computation of 
this expectation (such a formula may not exist for more complex dynamics of 
the underlying asset, or for exotic options), we will use a numerical approx- 
imation so as to illustrate our methodology on a realistic model example. 
The expectation is approached by the following Monte-Carlo procedure: 

r]{uo,K,d,r,^,R,So,T,K) = -^^^{Stj - K)+ 

i=i 

with an Euler-Maruyama approximation of iSt,i't)te[o,T] with timestep h > 
0: for j = 1, . . . , M and t = 1, r//i: 

So,j = So 

' Sthj = S(t-i)hj {l + Rh + ^u^t-i)h,jVh^Wl^>j 

^ vth,j = V{t-i)Kj (l + <0 - ^{t-i)Kj)h + C,/i7^njhjVh{rAWlj + Vl-r^^W^^j)) 

where /i > is the time discretization parameter, and {AW^j,A.W^j) are 
indepedent realizations of a standard Gaussian random variable. 
We fix Sq = 60 (the initial price of the asset), T = 0.25, K = 30, as well 
as the discretization parameters h = .001 and M = 10000. The uncertain 
parameters are X = (fg) and Z = (k, 9, r, ^, R), which are given the uniform 
distribution probabilities summarized in Table [H 

The coarse approximation uses a reduced number m of simulated trajectories 
to compute the empirical mean: 

-RT ™ 

Vcii^o, At, 9, R, t R, So,T, K) = V(Srj - K)+ 
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Note that for computing rjc, the same time discretization parameter h, as 
well as the same simulated Brownian increments AW'^''^ are kept, hence our 
approximation is "hierarchical" in the sense of Subsection 11.31 
We chose m = 5000, so that p = m/M = 1/2. 

2.2 Results and discussion 

We estimated ctc and o"e based on a sample of n = 100 realizations of each 
variable Y,Y',Yc and Y^. The estimates are: 

ac = .9017 ae = .4909. 

For comparison purposes, we also estimated 0"^,-^: 

a^n = -8491. 

We are interested in the (estimated) relative efficiency of the confidence 
intervals based on our variance-reduced estimator, as compared with those 
based on T'^, that is: 

ClassicalCost 

where ClassicalCost is the cost of the evaluations necessary to produce an 
asymptotic confidence interval of fixed length L using only the T'' estimator: 

ClassicalCost = 2—^{q(a)aTn)'^- 

As the denominator and the numerator of Eff are proportional to , the rel- 
ative efficiency is independent of the target length of the confidence interval 
L. 

In Figure [H we plot the estimated relative efficiency of our variance-reduced 
estimator, as function of the target risk level a. 

We see that, based on empirical estimations, our variance reduction enables 
an interesting reduction of the computational cost by more than 50% for 
a = 0.05, and this reduction is even more significative for small risk levels 
(up to 90% for a = 0.0001). 
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